Scale separation in granular packings: stress plateaus and fluctuations 
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It is demonstrated, by numerical simulations of a 2D assembly of polydisperse disks, that there 
exists a range (plateau) of coarse graining scales for which the stress tensor field in a granular solid 
is nearly resolution independent, thereby enabling an 'objective' definition of this field. Expectedly, 
it is not the mere size of the the system but the (related) magnitudes of the gradients that determine 
the widths of the plateaus. Ensemble averaging (even over 'small' ensembles) extends the widths of 
the plateaus to sub-particle scales. The fluctuations within the ensemble are studied as well. Both 
the response to homogeneous forcing and to an external compressive localized load (and gravity) 
are studied. Implications to small solid systems and constitutive relations are briefly discussed. 

PACS numbers: 45.70.Cc, 46.65.+g, 61.46.-w 



Continuum descriptions of matter comprise equations 
of motion for appropriate sets of macroscopic fields Q. 
It is helpful (though not required) when these fields or 
the constitutive relations do not depend on the averag- 
ing scale within a certain range of scales (larger than 
microscopic and smaller than the scales characterizing 
field gradients), see e.g., 0- These 'plateaus' of scales 
define the macroscopic fields and their fluxes in a "scale 
independent" way. The corresponding constitutive re- 
lations are local (expressed in terms of gradients of the 
fields) when scale separation exists. The Navier-Stokes 
equations and solid elasticity are typical examples. 

Recent studies render support to the notion of scale 
separation in nanoscale solids and granular matter, as 
their mechanics lends itself to description by elasticity 
for a range of loads IS H H S H E 113 ■ The short 
correlation length of the contact forces in static granular 
matter suggests that such plateaus should indeed 
exist. In this Letter, we study the scale dependence of 
the stress response in polydisperse granular packings. We 
show that even in small systems where stress gradients 
may be large, stress plateaus can be identified (though 
they may be quite narrow in some cases); expectedly, 
the widths of the plateaus are related to the gradients 
of the stress field. Ensemble averaging increases these 
widths. The scale dependence of the fluctuations within 
the ensemble is studied as well. 

The considered model is a two-dimensional (2D) rect- 
angular slab in the x — z plane, periodic in the x 
direction. It comprises 3600 polydisperse frictional 
disks, whose radii are uniformly distributed in the range 
(-Rmin, -Rmax = 2i?inin). The aspcct ratio is about 7. 
Gravity acts in the — z direction. It is prepared by se- 
quentially dropping the grains from rest at random hor- 



izontal positions above the top of the system. The floor 
comprises about 160 non-touching grains, whose centers 
are constrained to reside at z = 0, and whose radii are 
randomly chosen from the range (1.2i?i„in, i?max), which 
ensures that bulk particles cannot percolate through the 
floor. Similar systems have been studied experimentally 
both in two and three dimensions 0, 0, 0, • 

The force model used in the simulations is essentially 
that of Cundall and Strack 14]: overlapping disks are 
coupled by both normal and tangential springs, of respec- 
tive stiffnesses, fc„ = 1.5 • 10^ jj!^^^ and fct = 0.5fc„, where 
(to) is the mean particle mass and g the gravitational 
acceleration. A linear viscous damping force (dashpot) 
acts in the normal direction (parallel to the line con- 
necting the centers of the disks) , the damping coefficient 
chosen to correspond to critical damping. The tangen- 
tial forces are limited by the Coulomb condition, with a 
coefficient of (both static and dynamic) friction, /i = 0.5. 
The simulation is run until it reaches a numerically static 
state |l5|. To this system we apply either a uniform ex- 
ternal load (at its top) or a vertical compressive force 
acting on one particle at the top of the system. In the 
latter case the displacement gradients are large near the 
point of application of the force and decay with distance 
from this point; the effects of the local gradients on the 
widths of the plateaus can thus be studied in the same 
system. The external load is linearly increased from zero 
to its final value, Fa, 'yci a, time comparable to the typical 
relaxation time to static equilibrium. The load is subse- 
quently kept fixed, and the system is relaxed to a new 
numerically static state. The system's response is linear 
in the magnitude of the force for loads not exceeding a 
few times (to) g, as in Ref. 0. The response obeys su- 
perposition when two external forces are applied, and is 
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reversible to the slow removal of the force |l5| . Below we 
specialize to the range of loads for which the response is 
linear. Findings for single realizations as well as 'ensem- 
ble averages' over different realizations of the disorder 
are presented. The results reported here pertain to the 
normal stress response of the system (at the floor), i.e., 
the difference between the floor stress for the loaded sys- 
tem and the (near-uniform) floor stress of the unloaded 
system. 

The stress field, tTQ,^(r), at point r, is given (without 
the kinetic stress term, which vanishes in the static limit) 
by the following exact expression, which is fully compat- 
ible with the general equations of continuum mechanics 
(for both static and time dependent states) |3. ll^lrRll8j |: 
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where i,j are particle labels, a, (3 represent Cartesian 



components, 



Tj , where 



is the center of mass 



of particle i, and is the force exerted by particle j on 
particle i. The coarse graining (CG) function, (f>(R), is 
a positive semidefinite normalized function, with a single 
maximum at R = 0, and width w (the CG scale). The 
sign convention here is that compressive stress is positive. 

The validity of the equations of continuum mechan- 
ics (unlike the constitutive relations) is not resolution 
limited [with Eq. ^ defining the stress field], see e.g., 
1^ 01 . The actual values of the stress depend on the 
choice of cj) and w. For 'large' CG scales the dependence 
of the stress on the choice of is weak for 'reasonable' CG 
functions, such as [HEi: (/-(R) = l/iTrw'^)H{w - |R|), 
where H is the Heaviside function, and the Gaussian, 
0(r) = -i^e-t!''!/"')' (in 2D). Gaussian CG functions 
yield smoother stress fields than the Heaviside function, 
for obvious reasons [13. The stress field of Eq. ((U corre- 
sponds to the Born-Huang jl^ (or Irving-Kirkwood [2ll |) 
formula in the limit of large CG scales. 

Experimentally UQ, the normal stress at the floor, 
azz, is given by: a^^^ix) = -\Y.ifizH{^ - \x - Xi\), 
where L is the CG length (gauge area in experiments), 
fiz is the force acting on the floor particle i and Xi is 
its position. An identical expression can be obtained 
from Eq. by substituting the anisotropic CG func- 
tion (/.(r) = ii/(-| -\x- x^\)5{z). 

As a theoretical basis for defining ensembles for gran- 
ular solids is lacking, one usually averages (with equal 
weights) the desired entities over randomly generated 
realizations, subject to some constraints, such as the 
construction protocol of the system. It is not a-priori 
clear whether these averages indeed yield typical val- 
ues. In order to study the effect of ensemble averaging, 
we prepared 10 different samples; this rather small en- 
semble was extended as follows. To each of the realiza- 
tions we applied an external load Fq = 1 (m) g to differ- 
ent particles at the top of the slab. We verified that the 



corresponding fluctuations were nearly statistically inde- 
pendent when these positions were separated by distances 
exceeding about half the height of the system. The CG 
response profiles were therefore ensemble averaged over 
ensembles of effective sizes, N^, up to = 110. The 
x-coordinate of the load defines a; = 0. 

The inset of Fig. presents the response of the system 
to uniform forcing at its top both for the case of vertical 
forces (uniaxial vertical stress) and oblique ones (a com- 
bination of vertical stress and horizontal shear). Forces 
of magnitude 7 (m) g/rit, were applied to each of the 
rit ^ 150 particles whose centers resided at z > 20.5 {d), 
where (d) denotes the average diameter of a particle. The 
response is quite flat for w > (d) even for a single real- 
ization, and flatter for an average over five realizations. 
In contrast, in the inhomogeneous case (with a localized 
force at x = 0), we obtain narrower plateaus. Fig. de- 
picts the response at a; = in this case. Even for a single 
realization, one observes a (narrow) plateau. The devi- 
ations from the plateau at large values of w are due to 
the large scale (macroscopic) spatial dependence of the 
average stress. For the case studied here, the mean re- 
sponse function (for w of a few particle diameters and 
A^e = 110; see Fig.0J) can be well fitted by a Gaussian of 
half- width, W ~ 17.5 (d) (except at the tail), hence (ig- 
noring small scale fluctuations) the convolution involved 
in the coarse graining process yields, for w >> d, a (ap- 
proximate) Gaussian of half-width \/W^ + vP, in agree- 
ment with the results obtained by direct coarse-graining. 

The presence of a plateau suggests that the stress is 'lo- 
cally homogeneous'; however, since the forcing is macro- 
scopically inhomogeneous, the width of the plateau de- 
pends on position, as shown in Fig. [5] In order to study 
this dependence in more detail, we deflne the plateau 
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FIG. 1: Mean response to a force at x = 0, vs. CG width, w, 
for different ensemble sizes, A^e- The CG function is Gaussian, 
and the unit of a^^ is Fo/l, where I « 160 (d) is the width of 
the system; h is the height of the system. Inset: the response 
to homogeneous forcing at 45° and 90° to the horizontal. 
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FIG. 2: Mean response at different a:;-coordinates at the floor 
vs. CG width, w, for an ensemble size A^e — 110. The thick- 
ened parts correspond to the plateaus; • marks wo- 



width Aw as the largest connected range of w for which 
\[crzz{w) /azz{wo)] — 1| < e, for a given tolerance e (max- 
imizing Aw over < Wq < 5{d)). Fig. O presents the 
dependence of Aw on a: for e = 5% and Ng — 110. The 
plateau widths can be rather small (a few particle diam- 
eters) in parts of the system, as may be expected consid- 
ering the inhomogeneity of the response in this relatively 
small system. The dependence of Aw on x is related to 
the shape of the macroscopic response. One expects Aw 
to be smaller where the response is less homogeneous. 
In Fig. O we plot the first and second derivatives of the 
response (calculated for w ~ 4{d)). In general. Aw ap- 
pears to be anti-correlated with the second derivative of 
the response, rather than the first derivative. The reason 
seems to be that the average over a symmetric segment 
around x, of a profile which is locally linear (on aver- 
age) and whose fluctuations are essentially uncorrelated, 
yields basically the same (coarse-grained) value at x al- 
most irrespective of the width of the segment, w. It is 
easy to verify that the size of the plateau should be given 
e \<Jzz\'-, the peak in Fig.Ocorresponds 

vanishes and there the plateau 
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to a region where 
size depends on the next even (fourth) derivative of azz 

The dependence of Aw on for several values of x 
is shown in the inset of Fig. |3| As expected, ensemble 
averaging smoothes small scale fluctuations: while for 
iVe = 1 the plateaus start at it; > (d) (see also Fig. ^| , 
they probably extend down to w = d^s Ng ^ oo. Note 
that not only the forces fluctuate among the realizations, 
but also the particle positions. The plateau widths seem 
to practically saturate already for Ne « 50, suggesting 
that the fluctuations in the ensemble are not strong. 

The wide distribution of contact forces in granular 
materials 0, may suggest the existence of strong 



FIG. 3: The width of the plateau. Aw, vs. the a;-coordinate, 
for an ensemble size — 110 and tolerance e — 5%, as 
well as the first and second derivatives of the response (for 
w — 4(d)). Inset: Aw vs. A'^^ at different i-coordinates. 



stress fluctuations. However, since the force correlations 
are rather short ranged 11], their fluctuations are well 
smoothed by spatial averaging. Denote the standard de- 
viation (in the ensemble) of the response by Act. The 
relative standard deviation, Aa zzl <y zzi as a function of 
the CG scale w, for different ensemble sizes, is presented 
in the inset of Fig. ^ With increasing N^, Aozzl^zz 
seems to saturate to a well-defined limit (for a given u;) . 
The relative stress fluctuations at fixed x seem to approx- 
imately satisfy: Aozzj^zz oc w~'^l'^ . While an explana- 
tion of this possible "scaling" is still lacking, we believe 
it is quite surprising that the fiuctuations at small values 
of w (clearly related to the force and particle position 
fiuctuations) share the same (approximate) scaling with 
those at large values of w, which sample large-scale fluc- 
tuations and the spatial variation of the stress. 

The dependence of the stress fleld (for w within the 
plateau) and its fluctuations on the absolute horizon- 
tal distance from the load, \x\, is presented in Fig. 0] 
The fluctuations (Act) seem to decay nearly exponentially 
with x, at least not too far from a; = 0. The decay length 
depends on w. However, the mean stress decays faster 
than exponentially at large x (basically as a Gaussian), 
in conformity with the linear elastic solution for this case; 
the sharp drop near |x| = 50 (d) is due to the fact that 
the response becomes negative over a small interval, as 
expected from linear elasticity 0, 0| . A similar study 
of a Lennard- Jones glass in the linear response regime, 
with a localized force applied in the interior (rather than 
the boundary) is presented in Fig. 7 of for 2D and 
Fig. 8 of ^3 fo'^ 3D systems (the dependence on w is not 
discussed in 0jE|)- In our case the mean stress (calcu- 
lated at the boundary) decays faster than exponentially. 
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while in [7], [lOj the stress is calculated in the bulk, and 
decays algebraically (both results are consistent with lin- 
ear elasticity). This renders the relative fluctuations, in 
our case, non-negligible at large distances, as opposed to 
the decay with distance found in 0, 0| ■ Our findings 
suggest that this decay only characterizes intermediate 
distances (or possibly a wall effect), whereas asymptoti- 
cally the relative standard deviation is finite (see Fig.0J. 

In summary, our findings indicate that the magnitudes 
of the local gradients of the macroscopic fields (stress in 
the above case) determine the widths of the plateaus. 
Although large gradients are expected in nonuniformly 
forced small systems, the size of the system is not the 
main factor that limits these widths. The mere existence 
of the plateaus (which can be as small as 3-5 particle di- 
ameters) suggests that continuum theories may be valid 
for granular and mesoscopic solid systems, but one may 
need to go beyond simple linear descriptions. The satura- 
tion of the results for small ensembles (e.g., 40-50 realiza- 
tions) due to the short range of force correlations suggests 
that appropriate constitutive relations can be derived for 
such systems. Although this Letter is restricted to the 
linear response regime, we believe that the plateaus will 
continue to exist even near fluidization (they do exist in 
the fluid regime), hence, the above approach should be 
relevant to a rather large range of loads. 
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FIG. 4: Mean response and its standard deviation vs. ]a:| 
at the floor, for several values of the CG width w, for an 
ensemble size A^e — 110. Inset; Relative standard deviation 
of the response at x = 0, vs. w, for different A^e. 
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